function [result] = speccalc(dfile, v, time1, time2, a3, thresh)
% compute power spectrum of v data from t1 to t2
% does a quadratic subtraction to remove mean trends and low freq drift
% from data first
% 1-2 Mar 1999. P. Manis
% basic results of analysis are stored in ther result structure
% as:
% result.filename = 'filename'
% result.records = 'frec - lrec'
% result.f0 = peak osc. freq from spectral analysis
% result.max = psd at f0
% result.lat = latency shift maximum (i.e., jump)
% result.isi = interspike interval for pre-latency records
%
%
disp(sprintf('speccalc with: t1= %6.1f  t2=%6.1f  a3 = %d  thr = %6.1f', ...
   time1, time2, a3, thresh))
t = make_time(dfile);
%[v, t] = datac('getv'); % get the data first
%dfile = datac('getdfile');
result.filename = dfile.filename;
result.records = sprintf('%d - %d', min(dfile.record), max(dfile.record));
rate = dfile.rate*dfile.nr_channel/1000;
fsamp = 1000/rate; % get sampling frequency
fco = 2000;		% cutoff frequency in Hz
wco = fco/(fsamp/2) % wco of 1 is for half of the sample rate, so set it like this...
%[b,a] = elliptic(4, 0.3, 40, wco, 'low'); % compute elliptic filter coefficients for low-pass, 4th order, pretty flat
% [b,a] = chebys(4, 0.3, wco); % chebyshev filter
%[b,a] = butterworth(4, wco, 'low'); % butterworth filter
if(wco < 1) % if wco is > 1 then this is not a filter!
   [b, a] = fir_win(8, wco); % fir type filter... seems to work best, with highest order min distortion of dv/dt...
   vsmo = DigitalFilt(b, a, v')'; % filter all the traces...
else
   vsmo = v;
end

v = vsmo;
t1=min(find(t>time1));
t2=min(find(t>time2));
tdel=round(10/rate);
spike_thresh = thresh;
[first_spike, first_isi, nr_spikes, spike_train] = find_spikes(v, rate, spike_thresh);
leftwin = 16;
fswl = round((first_spike-leftwin)/rate);
fswr = round((first_spike+1)/rate);
slpos = round(((first_spike-2)/rate)); % spike latency position
t0a=logspace(0, 2, 26);
t0=round((t0a+20)/rate);
nr = size(v);
nfft =2048;
Pvgrand = zeros((nfft/2)+1, 1);
h = findobj('Tag', 'SpecPlot'); % check for pre-existing window
if(isempty(h)) % if none, make one
   h = figure('Tag', 'SpecPlot', 'Name', 'Spectral Analysis of Current Clamp', 'NumberTitle', 'off');
end
figure(h); % otherwise, select it
clf; % always clear the window...
%f = (1000/rate)*(0:(nfft/2))/nfft;
subplot(2,2,1);
title(sprintf('Spectral Analysis for %s Records %d-%d', dfile.filename, min(dfile.record), max(dfile.record)));
ylabel('mV');
xlabel('ms');
hold on;
Pv=zeros(nr(1), nfft);
ngrand = 0;
novlap = round(2/rate);
windlen = round(50/rate);

for i = 1:nr(1)
   if(slpos(i) < t2 & slpos(i) > 0)
      ts = slpos(i);
   else
      ts = t2;
   end
   if(a3 == 2)
      j = floor((i-1)/4)+1;
      t1l = t1+t0(j);
   else
      t1l = t1+tdel;
   end
   if(ts-t1l > (10/rate))
      tn = t(t1l:ts);
      vn = v(i,t1l:ts);
      p = polyfit(tn, vn,2);
      vl = vn - (p(1)*tn.^2 +p(2)*tn+ p(3));
      %     p = polyfit(tn, vn,1);
      %     vl = vn - (p(1)*tn +p(2));
      %     vl = vn;
      plot(tn, vl);
      [Pv, f] = spect(vl, nfft, novlap, windlen, 1000/rate);
      %vf = fft(vl, nfft);
      %      Pv(i,:) = vf.*conj(vf)/nfft;
      size(Pv);
      Pvgrand = Pv(:,1)+Pvgrand;
      ngrand = ngrand + 1;
   end
end
axis([150 200 -3 3]);
hold off;
Pvgrand = Pvgrand/ngrand;
   subplot(2,2,2);
   loglog(f, sqrt(Pvgrand(1:((nfft/2)+1))), '-k');
   axis([1, 10000, 0.01, 10]); % fix axis range
   [m,k]=max(Pvgrand(1:((nfft/2)+1)));
   ylabel('mV/Hz^{1/2}');
   xlabel('F (Hz)');
   text(f(k), Pvgrand(k), sprintf('\\uparrow\nFo = %8.1f Hz\nPer = %8.1f ms\nMax = %8.2f mV\^2/Hz', ...
      f(k), 1000/f(k), sqrt(Pvgrand(k))),...
      'FontSize', 6, 'VerticalAlignment', 'Top', 'HorizontalAlignment', 'Center');
   mx = max(Pvgrand(1:(nfft/2)+1))*2;
   %axis([10 10^3 10^-2 10^2]);
   subplot(2,2,3);
   hold on;

% store these results
result.f0 = f(k);
result.max = sqrt(Pvgrand(k));

v0=v(1,fswl(1));
%axis([-leftwin 2 v0-5 v0+10]);
for i=1:nr(1)
   %   subplot(2,2,3);
   hold on;
   if(slpos(i) < t2 & slpos(i) > 0)
      ts = slpos(i);
   else
      ts = t2;
   end
   if(a3 == 2)
      j = floor((i-1)/4)+1;
      t1l = t1+t0(j);
   else
      t1l = t1+tdel;
   end
   if(ts - t1l > 50) 
      tn = t(t1l:ts);
      vn = v(i,t1l:ts);
      p = polyfit(tn, vn,2);
      vl = vn - (p(1)*tn.^2 +p(2)*tn+ p(3));
      dvdt=diff(vl)./diff(tn);
      vx=vl(1:length(vl)-1);
      %      plot(vx, dvdt);
      %plot(t(t1l:ts-1), vx);
      %   plot(t(fswl(i): fswr(i))-(fswl(i)*rate)-leftwin, v(i,fswl(i):fswr(i)));
   end
end
ylabel('mV');
xlabel('t (ms)');
hold off;
subplot(2,2,3); % plot fsl vs ipi duration

fsl = zeros(1, nr(1)/4);
for i=1:4:nr(1)
   if(a3 == 2)
      j = floor((i-1)/4)+1;
      t1l = t1+t0(j);
   else
      t1l = t1+tdel;
   end
   % compute mean of fsl's in the block
   fsl(j) = mean(first_spike(i:i+3))-(t1-1)*rate-t0a(j);
   fsls(j)=std(first_spike(i:i+3));
   %   loglog(t0a(j), first_spike(i)-t1*rate-t0a(j), 'rs');
end
hlat = loglog(t0a, fsl);
set(hlat, 'Color', 'k', 'Marker', 's', 'Markersize', 3, 'MarkerFaceColor', 'k');
axis([1, 100, 1, 100]); % fix axis range
xlabel('Prepulse Duration (ms)');
ylabel('Latency (ms)');


subplot(2,2,4);
fsld = diff(fsl);
t0d = diff(t0a);
dldt = fsld./t0d;
[m, k] = max(dldt);
hdl = semilogx(t0a(2:length(t0a)), fsld);
axis([1, 100, -5, 15]); % fix axis range
text(t0a(k), dldt(k), sprintf('\\downarrow dl/dt = %8.1f t = %8.2f ms', m, t0a(k)),...
   'FontSize', 6, 'VerticalAlignment', 'Bottom', 'HorizontalAlignment', 'left');
xlabel('Prepulse Dur(ms)');
ylabel('dlatency/dt');

% store more result
result.lat = fsld(k); % derivative
result.jump=m; % largest jump
kk = (k-1)*4;
if(kk>0)
   result.isi = mean(first_isi(kk:kk+3)); % mean of isi 
else
   result.isi = NaN; % no isi...
end

return
